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Abstract 

Due to the long-term fluid-solid interactions in waterflooding, the tremendous variation of oil reservoir formation 
parameters will lead to the widespread evolution of preferential flow paths, thereby preventing the further enhancement of 
recovery efficiency because of unstable fingering and premature breakthrough. To improve oil recovery, the 
characterization of preferential flow paths is essential and imperative. In efforts that have been previously documented, 
fluid flow characteristics within preferential paths are assumed to obey Darcy's equation. However, the occurrence of non- 
Darcy flow behavior has been increasingly suggested. To examine this conjecture, the Forchheimer number with the inertial 
coefficient estimated from different empirical formulas is applied as the criterion. Considering a 10% non-Darcy effect, the 
fluid flow in a preferential path may do experience non-Darcy behavior. With the objective of characterizing the preferential 
path with non-Darcy flow, a hybrid analytical/numerical model has been developed to investigate the pressure transient 
response, which dynamically couples a numerical model describing the non-Darcy effect of a preferential flow path with an 
analytical reservoir model. The characteristics of the pressure transient behavior and the sensitivities of corresponding 
parameters have also been discussed. In addition, an interpretation approach for pressure transient testing is also proposed, 
in which the Gravitational Search Algorithm is employed as a non-linear regression technology to match measured pressure 
with this hybrid model. Examples of applications from different oilfields are also presented to illustrate this method. This 
cost-effective approach provides more accurate characterization of a preferential flow path with non-Darcy flow, which will 
lay a solid foundation for the design and operation of conformance control treatments, as well as several other Enhanced Oil 
Recovery projects. 
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Introduction 

During long-term waterflooding, the pore structure of a poorly 
consolidated sandstone formation will change significantly because 
of fluid-solid interactions, such as sand production [1] and clay 
erosion [2] . Therefore, the preferential flow path, which has also 
been termed the high-permeability streak or thief zone, is 
commonly formed in reservoirs because sand or clay particles 
have been eroded by fluid-induced force [3]. Existence of 
preferential flow paths will lead to inefficient water injection, 
unstable fingering and premature breakthrough of polymer or 
steam, all of which will prevent the recovery efficiency from being 
further enhanced. Accordingly, the identification and character- 
ization of preferential flow paths have become imperative for the 
Enhanced Oil Recovery process in reservoir development [4]. 

Approaches that can be applied to characterize the preferential 
flow path involve interwell surveillance [5-8] , well logging [9, 1 0] , 
reservoir engineering method [11-14] and well testing [15—18]. 
Darcy's law, which is the basis and principle of almost all of these 
methods, is always employed to depict fluid flow behavior in the 
preferential path. However, the influence of non-Darcy effect, 
especially for the high-rate wells with an ultra-high permeability 



streak, has been increasingly suggested for fluid flow within this 
path [18-20]. 

Our objectives in this work can be summarized as follows: (1) to 
examine whether non-Darcy flow behavior can occur in a 
preferential path; (2) to provide a hybrid analytical/numerical 
model to investigate the pressure transient response, which couples 
the non-Darcy effect of the preferential path with the Darcy flow 
of the reservoir; and (3) to present an available approach to 
characterize the preferential path with non-Darcy behavior using a 
pressure transient test. 

Theoretical Consideration 

Accurately describing the fluid flow phenomenon within the 
preferential path is essential for precise characterization and 
successful treatment project design. Flow through a preferential 
path is usually assumed to obey Darcy's equation, i.e., the pressure 
gradient is linearly proportional to the superficial velocity. 
Nevertheless, it's increasingly recognized that the non-Darcy 
effect should not be ignored because a preferential path always has 
an extremely strong transport capacity [4,18-20]. For example, 
tests from several poorly consolidated reservoirs in Chinese 



PLOS ONE | www.plosone.org 



1 



December 2013 | Volume 8 | Issue 12 | e83536 



Modeling of Preferential Path with Non-Darcy Flow 



Table 1. A brief summary of criteria for non-Darcy flow in porous media. 





Parameter 


Source 


Critical Value 


Investigation Method 


Re= P -^ 

V- 


Chilton and Colburn (1931) [35] 


40-80 


Experiments on packed particles 




Fancher and Lewis (1933) [36] 


10-1000 for unconsolidated, 
0.4-3 for loosely consolidated 


Crude oil, water and air through unconsolidated sands, lead shot, 
and consolidated sandstones 




Tek (1957) [37] 


1.0 


Air, water flow through Woodbine, Wilcox, Warren and 3 rd 
Venango sands 




Bear (1972) [38] 


3-10 


Review and analysis 




Scheidegger (1974) [39] 


0.1-75 


Review and analysis 




Dybbs and Edwards (1984) [40] 


1-10 


Experiments in fixed beds of arranged spheres and cylinders 




Blick and Civan (1988) [41] 


100 


Simulation of capillary-orifice model 




Du Plessis and 
Masliyah (1988) [42] 


3-17 


Representative unit cell simulation 


P 


Green and Duwez (1951) [43] 


0.1-0.2 


N 2 flow experiments through four different porous metal samples 




Ma and Ruth (1993) [44] 


0.005-0.02 


Diverging-converging model simulation 




Andrade et al. (1999) [25] 


0.01-0.1 


Simulation in disordered porous media 




Zeng and Grigg (2006) [23] 


0.11 


Review and theoretical analysis 




Chukwudozie 
et al. (2012) [45] 


0.02-0.08 


Lattice Boltzmann simulation of 3D realistic image-based porous 
media 




Ergun (1952) [46] 


3-10 


Gas flow experiments through packed particles 


ft l-<ji 








Re= P ^- 


Ma and Ruth (1993) [44] 


3-10 


Diverging-converging model simulation 


Re=^l 
f 


Thauvin and Mohanty 
(1998) [47] 


0.11 


Simulation of a pore network model 


Rc - 4pTv 


Comiti et al. (2000) [48] 


4 


Theoretical prediction and experimental data published previously 


p fia,j(l-(j 


) 






p\/kv 

Ret = 

/' 


Newman and 
Yin (2011) [49] 


0.1-0.3 for cubic flow, 1-3 for 
Forchheimer Equation 


Lattice Boltzmann simulation of stochastically generated 2D 
porous media 


In this table, a vd is the dynamic specific surface area, m -1 ; d t is the throat diameter, m; r is the radius of pore throat, m; 7" is the hydraulic tortuosity, dimensionless; u is 
the fluid intrinsic velocity, m/s; <I> is the porosity, dimensionless. 
doi:1 0.1 371 /journal.pone.0083536.t001 



oilfields indicate that the effective permeability of interwell 
formations can reach up to 10-150 (im 2 , and tracer advance 
speed can be as high as 10-100 m/h, which is absolutely 
incredible for Darcy flow [21,22]. Therefore, in this section, using 
the criteria proposed by Zeng and Grigg [23], non-Darcy flow in 
the preferential path is identified with tracer test data from the 
Gudong Oilfield. 

Darcy's law is widely applied to characterize fluid flow in porous 
media. However, a large number of experimental observations 
[24] and numerical simulations [25,26] have demonstrated that 
Darcy's law should be restricted to viscous flow conditions or, 
more precisely, to low values of the Reynolds number (Re), which is 
given by Eq. (1). 



Re- 



(1) 



In this formula, p is the fluid density, kg/m 3 ; v is fluid superficial 
velocity, m/s; Dp is the grain diameter, m; and fi is the fluid 
viscosity, Pa-s. 

As the Reynolds number increases, due to the contribution of 
inertial forces to laminar flow through the void space, the 
distribution of streamlines along the direction orthogonal to the 



main flux becomes more homogeneous, which causes flow 
deviations from Darcy's law [25]. The deviation from the classical 
linear relationship is defined as non-Darcy flow, which can be 
quantified macroscopically by an empirical formula named the 
Forchheimer equation [27], 



dp 
dx 



(2) 



where p is the fluid pressure, Pa; x is the space coordinate in the 
fluid flow direction, m; k is the permeability of the porous media, 
m 2 ; ft is an experimentally derived parameter called the non- 
Darcy flow coefficient or inertial coefficient. Non-Darcy flow 
behavior has been recognized as a common and important 
phenomenon in the development of hydrocarbon reservoirs and 
geothermal resources [28], especially for the fractured wells [29] 
and the near-well region of high-rate gas wells [30]. The 
significant influences of nonlinear flow on the well productivity 
and reservoir performance have been analyzed and confirmed 
[31,32]. 

Inconsistent with the sudden passage from laminar to turbulent 
flow in conduits and channels, where there is a critical condition 
separating both regimes, the transition from linear to nonlinear 
flow behavior is more likely to be gradual for the transport in 
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Water Cut 



Water Cut 



Figure 1. Formation parameter variation for the Gudong Oilfield during long-term waterflooding. 

doi:1 0.1 371 /journal.pone.0083536.g001 



porous media [25,33,34]. To identify the onset of non-Darcy flow, 
several criteria have been proposed in previous work [23,25,35- 
49], which can be summarized in Table 1. Table 1 shows that 
there are roughly three types of criteria for the determination of 
inertial flow beginning through porous media: Type 1 represented 
by the Reynolds number (Re), and the critical value varies from 0. 1 
to 100; Type 2 by the Forchheimer number (Fo) ranging from 
0.005 to 0.2; Type 3 by various types of modified Reynolds 
numbers (MRes), and the critical value varies with different 
definitions. Since the Reynolds number has its root from a similar 
criterion for turbulent flow in pipes, both Type 1 and Type 3 have 
been employed primarily for packed particles and unconsolidated 
or poorly consolidated sands, for which a characteristic length is 
available. Unfortunately, without a well-defined characteristic 
length, these criteria cannot be employed. However, as recom- 
mended by Zeng and Grigg [23], the Forchheimer number (Fo), 
which is defined by 



Fo- 



kfipv flpv 2 



(3) 



is of explicit definition, rational physical meaning, and general 
applicability [50]. Eq. (3) illustrates that the Forchheimer number 
represents the ratio of the pressure gradient consumed by 



fluid-solid interaction ppv 2 to that by viscous resistance, (xv/k. 
For any types of porous materials, as long as the permeability and 
non-Darcy coefficient are obtained either by experimental 
techniques or by empirical models, the Forchheimer number 
can be available. Thus, with interwell tracer test data from the 
Gudong Oilfield, the Forchheimer number is employed as the 
criterion to determine whether non-Darcy flow will occur in a 
preferential path. 

Fig. 1 shows changes of formation parameters from core 
analysis of the Gudong Oilfield, which has been developed by 
waterflooding for 26 years. From Fig. 1, it can be observed that 
along with the increment in the water cut, the permeability, 
porosity, and pore/throat radius of all the layers are experiencing 
a dramatic increase, whereas the content of shale is decreasing. 
Taking the formation permeability of Layer 4 as an example, at the 
initial stage (water cut = 10%), its permeability is only 2000 x 10 3 
um 2 . However, when the water cut increases to 90%, the 
permeability has become 4928x10 J u.m 2 , which is almost 2.5 
times as large as the initial value. Analysis of injector/producer 
production performance also demonstrates that large numbers of 
preferential flow paths (confirmed by the interwell tracer test) have 
evolved in this oilfield. Recent statistics from 56 well pairs in this 
field reveal that in 82% of the tracer tests, the tracer advance speed 
can reach as much as 15-82.5 m/h, i.e., the tracer can be detected 
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Table 2. Empirical models for non-Darcy coefficient estimation. 



No. 



Source 



Empirical model 



Investigation Method 



Ergun (1952) [46] 

Janicek and Katz (1955) [52] 
Geertsma (1974) [53] 
MacDonald et al. (1979) [54] 

Pascal et al. (1980) [55] 

Jones (1987) [24] 

Coles and Hartman (1998) [56] 

Coles and Hartman (1998) [56] 

Li et al. (2001) [57] 

Friedel and Voigt (2006) [58] 

Friedel and Voigt (2006) [58] 



100 /l.8xl0» 



1.82 x 10 10 
1.59 x 10 5 

100 /245x 10 8 



4.8 x 10 12 



2.018 x 10" 

¥^ 

3.51 x lO 1 ^"- 449 

¥™ 

8.17 xlO"/ 537 

£1.79 

1.15x10' 



k<f> 
4.1 x 10 

1.1 x 10 



11 



C0 2 , N 2 , CH 4 and H 2 flow through various sizes of spheres, sands, and pulverized 
coke 

Flow through sandstone, limestone and dolomite 

Both consolidated and unconsolidated sandstone, limestone and dolomite 

Experimental data from previous work, including spherical glass beads, cylindrical 
fiber beds and consolidated media 

Multirate field test of low permeability hydraulically fractured wells 

Experiments of vuggy limestone, crystalline limestone and fine-grained sandstone 

Sandstone and limestone samples for flow testing using the same porosity method 

Sandstone and limestone samples for flow testing using the simultaneous 
equations 

Numerical simulation of N 2 flowing at various rates through wafer-shaped Berea 
sandstone 

Experimental data from previous studies 

Experimental data from fractures with a variety of proppants 



In formulas cited above, k appears in 10 3 jim 2 . 
doi:1 0.1 371 /journal.pone.0083536.t002 



from the production liquid of observation wells several hours after 
it was injected [22]. 

To determine whether non-Darcy behavior has emerged in a 
preferential flow path using the Forchheimer number, estimating 
the inertial coefficient is essential. Due to the unavailability of 
experimental data, an empirical model is applied to calculate the 
non-Darcy coefficient, which assumes a direct correlation of this 
parameter to fundamental properties of the porous media, such as 
permeability and porosity. Although tortuosity is also a key factor 
influencing non-Darcy behavior and is included in some 
correlation formulas [51], the tortuosity of the preferential flow 
path is unavailable; therefore, empirical models with tortuosity are 
not used here. 

Table 2 provides a sample of correlations for non-Darcy 
coefficients from previous studies [24,46,52-58]. By substituting 
data from the block of interest (£=5u.m 2 , <J> =35.5%), inertial 
coefficients are estimated for each model and shown in Table 3. 
The values of non-Darcy coefficients obtained with different 
formulas are quite distinct from one another, which makes it 
difficult to provide definite interpretation. According to the 
interwell tracer test data from the Gudong Oilfield, with each 
listed in Table 3, the Forchheimer number is calculated for 
different tracer advance speeds ranging from 15 m/h to 82.5 m/h 
(Fig. 2). The red line of critical Forchheimer number 0.11 is also 



plotted in this figure, corresponding to a 10% non-Darcy effect 
[23,59]. Comparison of these plots in Fig. 2 reveals that for almost 
all of the inertial coefficients estimated with these formulas (except 
No. 8 proposed by Coles and Hartman (1998)), fluid flow above a 
critical velocity (v c ) within a preferential path will experience non- 
Darcy behavior. Referring to the inertial coefficients in Table 3, 
larger \i values are observed to lead to smaller critical velocities. 
This picture has therefore demonstrated that if the fluid velocity 
equals or exceeds the critical value, non-Darcy flow will occur in 
the preferential flow path, consistent with experimental research of 
Liu et al. [20]. Description technology for the preferential path 
with Darcy flow can be found in previous studies [5—18]. 
However, no documented effort to take the non-Darcy flow effect 
into account in the characterization of a preferential flow path has 
been observed. 

Methods 

Physical Model and Assumptions 

Preferential flow path is commonly believed to be a tubular 
porous medium with extremely high permeability and large pore/ 
throat radius, in which fluid flow is dominant for the whole 
reservoir [4,13]. Fig. 3 presents the diagrammatic representation 
of a vertical well intersected by a preferential flow path. Although 



Table 3. Non-Darcy coefficient estimated by different empirical formulas. 



Formula No. 



10 



11 



Value of II, 10 5 m ' 2.8367 



doi:1 0.1 371 /journal.pone.0083536.t003 



PLOS ONE | www.plosone.org 



December 2013 | Volume 8 | Issue 12 | e83536 



Modeling of Preferential Path with Non-Darcy Flow 



10 



0.1 



50.01 



/c=5um 2 ,(#=0.355,/^1000kg/m 3 ,//=0.2mPa.; 
-■—No. J — •— No. 2 -^*^No.3 
~*-No.S — «— N0.6 — *-No.7 -»-No.8 
— «— No.9 -»-No.10-»-No.11 




Fo=0.11 



10 20 



30 40 50 60 
Fluid velocity, m/h 



70 



80 90 



Vertical Well 



Preferential Flow Path 



y 



Figure 2. Calculated Forchheimer number for tracer advance 
speed ranging from 15 m/h to 82.5 m/h and the comparison 
with the non-Darcy flow criterion. 

doi:1 0.1 371 /journal.pone.0083536.g002 

the mathematical model proposed in this work is completely 
general and can be extended to various combinations of boundary 
conditions, we will concentrate on an infinite slab reservoir with 
impermeable upper and lower strata. 

The compressibility of liquid in the formation can be attributed 
mainly to the content of gas, i.e., free gas and solution gas [60]. 
Predetermined by the geological conditions, less dissolved and free 
gas is contained in the liquid of poorly consolidated or 
unconsolidated sandstone formations, where the preferential flow 
path is widely formed. Moreover, the preferential flow path always 
evolves after long-term waterflooding (approximately 10 to 
30 years). At that time, the gas has been almost totally exploited 
such that the content of gas is very small or approaches zero [61]. 
Thus the oil can be considered slightly compressible [16,18,62- 
65]. 

Compared with several decades of reservoir development, the 
pressure transient test period is very short (less than seven days), 
and the intrinsic topological structure of the reservoir formation 
hardly experiences tremendous variations. The shear therefore 
does not change during the test and the influence of the shear on 
the viscosity can be negligible [66] . Although the fluid viscosity is 
also affected by other physical parameters, such as temperature 
and pressure, these parameters can also be considered to be 
invariable during the short test. Hence, assuming that the fluid is 
of constant viscosity is reasonable [16,18,62-66]. 

To sum up, the following assumptions are made to facilitate this 
mathematical problem: 

(1) The reservoir is filled with a single-phase, isothermal, slightly 
compressible fluid, for which the compressibility and viscosity 
are assumed to be constant. 

(2) The formation is of two-dimensional anisotropy, and k, k v 
stand for the horizontal and vertical permeability, respective- 

iy- 

(3) The bottom (z = 0) and top (z = h) boundaries of the formation 
are impermeable, and the lateral boundaries do not manifest 
themselves during the time period of interest. 

(4) The preferential flow path extends horizontally in the x 
direction and is located at the elevation Zp from the lower 
boundary of the formation. 

(5) Due to its small volume, one-dimensional (axial) non-Darcy 
flow is assumed within the preferential flow path. 



Figure 3. Schematic diagram for a vertical well intersected by a 
preferential flow path. 

doi:10.1371/journal.pone.0083536.g003 

(6) At the intersection point of the preferential flow path with this 
well (x = 0), the production rate q is assumed to be constant, 
whereas at the other end {x = L), the flow rate is 0. 

(7) At different times and different locations, the flow rate from 
the reservoir to the preferential path is different. Fluid enters 
the preferential flow path along its axial direction at an 
unspecified rate per unit length, q k (x, t), and the flow rate at 
any point x within the preferential flow path is qjx, t), then 



qi, [x ,t\dx 



q d (x,t) -- 



Dimensionless Variables 

Dimensionless variables are defined as follows: 
(1) Dimensionless pressure drop: 

kh(pi -p) 



Pd ~- 



(2) Dimensionless time: 



1.842 x 10- V< 



t D - 



(3) Dimensionless distance: 



3.6 x Wkt 



x y z p 

XD -Lfr yD ~Lfr ZpD ~^ 



(4) Dimensionless length of the preferential flow path 



" 2/7 V k 

(5) Dimensionless radius of the preferential flow path: 



(4) 



(5) 



(6) 



(7) 



'pD '- 



h 2h 



1 1 
k\4 , fk v \i 

k v 



(8) 
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(6) Dimensionless flow rate: 



Id IhL 
IdD = — qhD, = — 
q q 



(7) Dimensionless conductivity of the preferential flow path: 



(9) 



CdD = 



K f nr p 

kh(L/2) 



(10) 



(8) Dimensionless flow rate constant: 



(9 



1 



DND) d z 



86400 nr 2 n 



(11) 



(9) Dimensionless storage coefficient: 



Cz>=- 



2Kij>C t h{L/2Y 



(12) 



where k and k v are the horizontal and vertical reservoir 
permeability, m ; h is the formation thickness, m; />,- is the 
initial reservoir pressure, Pa; q is the flow rate at the heel of the 
preferential path, m / d; qh is the influx rate per unit length, 
m3/(d-m); qd is the cumulative flux at a point in the 
preferential flow path, m3/ d; t is the production time, hours; 
<J> is the reservoir porosity, dimensionless; C t is the total 
compressibility, Pa -1 ; L is the length of the preferential flow 
path, m; Zp is the elevation of the preferential flow path from 
the lower boundary, m; rp is the radius of the preferential flow 
path, m; Tp is the quivalent radius of a preferential flow path in 
an anisotropic reservoir, m; kf is the permeability of the 
preferential flow path, m ; C is the storage constant, m 3 /Pa; x, 
y, z are the space coordinate in the x, y, ^-direction (m), 
respectively. 



Mathematical Model 

Mathematical model of non-Darcy flow in preferential 
path: Fluid flow in a preferential path with non-Darcy behavior 
can be represented by the Forchheimer equation: 



ax kf 



(13) 



where pjx, t) is the pressure in the preferential flow path at point * 
at time t, Pa. 

Fluid velocity in the preferential flow path at point x at time t is 



qd(x,t) 
86400ot-2 



Substituting Eq. (14) in Eq. (13), we have 



(14) 



Spd 

dx 



H q d {x,t) 



86400OT-2 



+ PP 



qd{x,t) 



86400W-2 



(15) 



Eq. (15) can be simplified by applying those aforementioned 
dimensionless variables 



dpdD 2n 271 2 

--^— = -^9dD+ 7 ^(qDND) d q dD (16) 

OXd <-dD ^dD 

Eq. (16), which describes the pressure drop per unit length along 
the preferential flow path, must be integrated to determine the 
dimensionless pressure at point xD. Integrating Eq. (16) for xD 
from 0 to xDj yields 



•*Dj 
0 



OPdD 

Sx D 



2k 

CdD , 

2k 



■*Dj 



0 



qmdx D + 



CdD 



l x ty 7 I 
{qDND) d qdD dx D 

Jo 



(17) 



Eq. (17) can be written as 



PwD —PdDj z 



2n [ XD J i 2n 

- q d D<lx D + — — (qDNDld 

bJo 



c d . 



CdD 



f J q 2 dD dx D 

Jo 



(18) 



where pw is the pressure at the heel of preferential flow path, Pa. 

Eq. (18) is an integral equation, which is very difficult to 
evaluate directly. Nevertheless, it can be calculated by discretizing 
the integral in space [67,68]. The preferential flow path is 
therefore divided into M equal segments (Fig. 4). The flux density 
is assumed to be constant in each segment, whereas the cumulative 
flow rate, q llD , varies from point to point along the preferential 
path, which can be approximated at point x' D as follows. 



qdD = qdDi + 



qdDi+l — qdDi 

Ax 



Ax 



(19) 



In this equation, Ax = 2AM is the size of each segment; i refers 
to segment i. 

From Eq. (18), the integrals can be evaluated separately, 



| (| f[x D Jdx D = j f[x D Jdx D + 

Jx Dj-\ + ~T 

t{Cl f ^ d * D Y 



(20) 



x Dj-^r 



f[x D )dx D 



The following integrals of the polynomials can be calculated 

Ax\ ./ 1 



f x Di+^2~ 



Jx d,-^t 



[ X Dj 
Dj 2 



X D -X D i+ ■ 



dXr 



-(Ax) 



x D -x Di +^)dx D =UAx 



(21) 



(22) 
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X D -0 X D j 



4 

< ► 


► 




• 

1 


• 

2 




< 


> 


• 




• 


< 


» 




• 

M-1 


• 



Figure 4. Diagram showing segments with cumulative flux distribution in the preferential flow path. 

doi:10.1371/journal.pone.0083536.g004 



A ( x D - x Di + — I dx B = - (Ax) 



^ / ' Ax\ , 1 /A x3 



By applying Eq. (19), one then obtains 
qdDdx' D = q dDi -Ax + 



qdDi+\—qdDi (Ax) 



Ax 



Ax 



(qdDi+l+qdDi) 



(23) 



(24) 



(25) 



X Dj 2 , > _ 2 Ax 

q dD dx D — q Wj ■ — + 
r „._Al Z 



„ qdDi+i-qdDi (Ax) 

2to ' Ax 'T + 



(28) 



/ gjgi+i -grfpA 2 (Ax) 
I Ax / 



"24~ "24 
1 



-AxqdDj-qdDj+\ + —Axq 2 dDj+l 



According to Eq. (20), the integrals in Eq. (18) can be expressed 



x Dj 



/-I 



q dD dx' D = ^2 (^qdDi 



Ax 



-qdDi+\ 



3 A 1 , 

+ - Ax- q dDJ + - Ax-q dDj+ 1 



(29) 



' X °J , i Ax qdDj+i-qdDj {Ax) 

qdD<lx D = qdDj ' -z- + 

x Dj 2 



Ax 



3 a 1 
-qdDj-Ax+ -qdDj+vAx 



v i Ax 

v Ax 
V B/ — 2" 



c l 2 dD dx D= C ldDi-^ X + 



2qd 



qdDi+i-qdDi {Ax) 



Ax 



f grn+i -qdDi \ 2 {Ax) _ Ax , 
I Ax J 3 " 3 £fe + 



Ax 



Ax 



„-qdDiqm+\ + -g-toi+i 



(26) 



(27) 



7-1 



Ax 
T 



2 , / v-i / ax , 

q dD ax D - I -z-qdm + —qdmqdDi+\ + —q~dDi+\ 



Ax 2 

T 



1 



+ 74 A - v ^»/ + &xqdDj-qdDj+ 1 



(30) 



24 
1 

24 



Now the differential equation of non-Darcy flow in the 
preferential flow path has been established by Eqs. (18), (29) and 
(30). 

Mathematical model of Darcy flow in the reservoir: If 

the preferential flow path is considered to be a line source, then the 
dimensionless pressure drop p dD at point x D on its surface can be 
obtained analytically with Green's function [69,70]. 



PdD{XD,r p D + ,to) 



0 J 



{xd — ocY 



4(to-T) 

{\+2^2exv[-n 2 n 2 L D {t D -T)} 



(31) 



[nn (z pD + r pD ) ] cos (mtZpB ) } 



dx 
t D -i) 



PLOS ONE | www.plosone.org 



7 



December 2013 | Volume 8 | Issue 12 | e83536 



Modeling of Preferential Path with Non-Darcy Flow 



The dimensionless pressure can be discretized on the surface of 
the preferential flow path using the same grid and time step as the 
non-Darcy flow equation to yield the following expression. 



•t D M i \ i 

PdD{xDj,r P D,tD)= qhDi (t D ) Si [xpjjp — t D j dt D (32) 



o 7=1 



where 



S,(x d ,t) -- 



~4~ 



I 1 2> 2 



,T/ U:l V " M 



{ 1 + 2 ^ exp [ -hVlJt] 



(33) 



cos[nn(z pD + r pD )] cos(nnz pD )} 

In this equation, erf is the error function. Eq. (32) can also be 
discretized in time as follows. 



IhDi \^D) ( x Dj> to t D \ dt D 



(34) 



Here, JV is the number of time step. If qux is assumed to be 
constant during the time interval (t Dk - t Dk - {), then Eq. (34) can be 
written as 



PdD{xDj,r p D,lD) =PdD{ X DjSpD,tDN) + 

J£, t'DN -<DN-1 

2, IhDi(tDN) 

i=i Jo 

where 



, ( 35 ) 
SAx D j,t D )dt D 



l TT^'v^ \c'DN-'Dk-l / , \ , C'DN-'Dk 

p, a ,(^j«.,M=i 1 i f «»»('*)y o s,(x DI ,t D )dt D -^ Si (x Bi j B )d; B 

(36) 

Coupling Relationship: By the continuity of pressure and flow 
rate on the surface of the preferential flow path, the p^ terms in 
Eqs. (18) and (35) are the same. The mathematical relationship 
between cumulative flow rate qjo and flux density q^o flowing from 
the reservoir to the preferential path at point x can be given by 



q d {x,t)=q- q h {x,t)dx 
Jo 

= q- | ° ^q h (x,t)dx D 

= q- \ \ D qqi,Ddx D 
1 Jo 

Eq. (37) can be expressed in dimensionless form as 



(37) 



IdD = 1 — ^ | QhodXr) 



Differentiating Eq. (38) yields 
dq d D 



, = - ~ qiiD 

dx D 2 



(38) 



(39) 



From the grid we used to discretize the preferential flow path, 
Eq. (39) can be approximated by 



qhD-- 



dq dD 2(q dDi — q dDi+i ) 



dxr 



Ax 



(40) 



Thus, substituting Eq. (40) into Eq. (35) yields 



PdDj —PdD { x DjS p dJdn) + 



2(qdDi — qdDi+i 



M 

E 

<dn-'dn-\ 



Ax 



(41) 



Jo 



i^XDjJo^jdtj) 



Solution Procedure 

Substituting Eq. (41) for the term p dI yj in Eq. (18), the following 
discretized form of non-Darcy flow equation within the preferen- 
tial flow path is obtained, 



' / . \ sr^^ldDi — qdDi+x) 
PwD-p dD \XDj,r pD ,t DN ) - ^ 

l 



'dn-'dn-i 



0 

271 
CdD 



>i(x D j,t' D ^dt' D 



^ — ■» / AX L\X \ D 

2^ I -yqdDi+ -yqdDi+l I + -Ax-q dDj + -Ax-q dDj+ 



Ax 



1% 

+ -7^(qDND) d 
<-dD 

1 



7-1 



^-^ / Ax 2 Ax Ax 2 

2_, y-y9dDi+ -^qtDiqdDM + -yq d Di+ 



1 



+ ^ Ax qdDj + g &xq dDr q dDj+ 1 + — Axq dDj+ , 



(42) 



Eq. (42) represents the dimensionless pressure loss between the 
vertical well intersection and point Xj^. Evaluation of Eq. (42) at 
the center of each segment xjyj brings about a set of M equations 
with M+2 unknown variables: p m r>, (jdDj (t= l,2,...,Af+l). Two 
more boundary conditions are therefore required. The first 
boundary condition can be given by the assumption that the flow 
rate at the tip far from this well is 0. 
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to(M+l)=0 



(43) 



Since the flow rate at the intersection point is equal to q, the 
second equation can be expressed as 



1 



(44) 



Eq. (42) can be solved in a forward manner in time because the 
evaluation of the p'jo term in Eq. (42) needs the knowledge of qjoj 
(i= 1,2,...,M+1) at old time steps to\,to2,--->tDjv-i- In addition, the 
set of equations is nonlinear, so the Newton-Raphson iterative 
scheme is employed at each time step. 

Moreover, the following equation can be used to incorporate 
the effect of storage Co [71] 



1-C, 



dp wD {t) 



ch 



dp„, D 



dt D 



to—i) 



dx (45) 



where p m is wellbore pressure without the storage effect, which is 
obtained from our model, Pa. The parameter p' w is wellbore 
pressure with the storage effect at the heel of the preferential flow 
path, Pa. 

Discretizing Eq. (45) yields the following formula [72] 



1 

] i C DPwD{'Dn- t Dn-\) 
t Dn-'Dn-\ 



9diPwb (tad- £ ?«+i)'P"J>( ( a»- f »+i) 

- qm - 2 P„D ( t Bn - i ) + 1 + C D ^g^gj 'Pwd {to, -'n„-i 



(46) 



where 



qm = 1 — Cj 



tDi — tDi-\ 



(47) 



Model Verification 

For the mathematical model, the appropriate choice of the 
segment number used to discretize the preferential flow path in 
space is essential to guarantee the accuracy and efficiency of the 
fluid flow simulation. Increasing the discrete element number M 
would improve the accuracy, whereas the size of the coefficient 
matrix needing to be solved also increases dramatically at each 
time step, leading to a long computation time. To provide a 
measure, we compare the pressure transient response of different 
segment numbers for the base-case model. Fig. 5 presents the 
effect of discrete element numbers on the dimensionless wellbore 
pressure. The result shows that if the number of segments is 
greater than or equal to 60, these pressure curves would 
completely coincide; however, if M is less than 60, there is a great 
difference from the other cases. Therefore, to ensure the accuracy 
and to reduce computational costs, a value of 60 segments is 
applied and recommended for the simulation. 

For the reason that no transient-pressure solution for a 
preferential flow path with non-Darcy behavior is available in 
previous literature, the model presented in this work is verified by 
examining the analogous case of a horizontal well, which couples 
wellbore hydraulic and reservoir fluid flow [70]. Because flux 
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Figure 5. Influence of discrete element numbers on the 
dimensionless pressure. 

doi:10.1371/journal.pone.0083536.g005 

distribution dictates both the flow behavior and the pressure 
performance, we will compare the flux distribution along the 
preferential flow path with the flux distribution of a horizontal 
well. Since the preferential path has an extremely high flow 
capacity, the physical model of the horizontal well is similar to the 
preferential flow path, except that fluid flowing through the 
wellbore is pipe flow instead of non-Darcy flow. From this 
perspective, characteristics of the influx profile should be similar 
for both the preferential flow path and the finite-conductivity 
horizontal well. 

Fig. 6A, showing the dimensionless influx distribution along the 
preferential flow path at different times, indicates that at early 
times (to — 10 4 ), most of the fluid enters the preferential flow 
path from the near vicinity of the wellbore, and flow across the tip 
far from the wellbore intersection is negligible. This observation is 
in agreement with theoretical analysis because pressure drop 
always starts from the near wellbore and transmits outwards, 
which can also be substantiated by the dimensionless flow rate 
distribution at different times (Fig. 6B). With the time increasing, 
formation far from the intersection begins to supply fluid 
gradually, so flux increases in most parts of the preferential flow 
path except the near intersection. Moreover, the largest increase of 
flux happens at the tip far from the vicinity of the wellbore. Fig. 6A 
also reveals that the flux distribution stabilizes at late times (to 
S10 3 ). At that time, the production rate comes mainly from the 
vicinity of the formation of the tips because of their larger drainage 
area. In addition, although the flux profile is a U-shaped curve, it's 
asymmetric with respect to the mid-point of the preferential flow 
path, which can be attributed to the influence of non-Darcy flow. 
Like the frictional pressure drop along a finite-conductivity 
horizontal well, inertial flow leads to a greater pressure drop at 
the near vicinity of the wellbore intersection, from which most of 
the fluid enters the preferential flow path. Therefore, we have 
attained good qualitative agreement with the influx profile of the 
finite-conductivity horizontal well, which was proposed by Ozkan 
et al. [70]. 

Results and Discussion 

In order to discuss the pressure transient response of a 
preferential flow path with non-Darcy behavior, we present the 
results obtained by the hybrid analytical/numerical model 
developed in this work. The influence of the following parameters 
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Figure 6. Dimensionless influx rate distribution at t D =0.0001, 0.2 and 1000 (A) and dimensionless flow rate distribution at t D 
= 10 4 , 10 3 , 0.02, 1.128, 14.38 and 1000 (B) along the preferential flow path (L D =2.5, r pD =0.05, C dD =5, z pD =0.5, q DND =5, C D =0). 
doi:10.1371/journal.pone.0083536.g006 



on the pressure performance of the preferential flow path is also 
discussed: length (L D ), radius (rj>n), conductivity {CdD), vertical 
location (^d), flow rate constant ((qDXDldi an d storage coefficient 
(Co). Then a computer-aided automatic matching algorithm is 
applied to characterize the preferential flow path using pressure 
transient analysis. 

Pressure Transient Behavior 

Fig. 7A shows the dimensionless pressure and derivative 
responses as a function of time for the preferential flow path 
without the storage effect. It can be observed that the general 
pressure transient behavior of this model can be divided into 
several flow periods. 

Stage 1: One-dimensional non-Darcy flow in a preferential 
path, which is analogous to early-time linear flow in a high- 
conductivity fractured well. However, because the preferential 
flow path, which evolves from long-term waterflooding, is always 
of extremely high conductivity, the flow is so rapid that the linear 
flow occurs too early to be of practical interest and is therefore not 
displayed in Fig. 7A. In practice, even a small storage coefficient 
(C B =5xl0" 4 ) may mask the characteristics of the early- time one- 



dimensional flow and the transition regimes, as indicated in 
Fig. 7B. 

Stage 2: At intermediate times, there is a vertical radial flow in 
the formation and one-dimensional non-Darcy flow within the 
preferential flow path. This is characterized by an approximately 
straight line and similar to the radial-linear flow in the horizontal 
well intercepted by a finite-conductivity vertical fracture (indicated 
by almost constant derivative response) [73] and the bilinear flow 
in the vertical well with a finite-conductivity vertical fracture 
(indicated by a 1/4-slope straight line) [74]. 

Stage 3: As time increases, a formation linear flow will be 
observed, which is indicated by a half-slope straight line for both 
the dimensionless pressure and derivative curves. This flow regime 
is akin to the linear flow in finite-conductivity vertical fracture 
communicating with a vertical well over its entire height. 
However, the gap between pressure P D and derivative dP D /dlnt D 
is not lg2 =0.3010 but approximately 0.225. Like the pressure 
response of the horizontal well, if the ratio of the preferential flow 
path length to the formation thickness is very large, this regime will 
last a very long time [75] . 

Stage 4: At late times, the system reaches a pseudo-radial flow 
period, and the pressure derivative curve becomes a horizontal 
line. However, the formation linear flow (Stage 3) always continues 




Figure 7. Pressure transient response without (A) and with (B) storage effect. 

doi:1 0.1 371 /journal.pone.0083536.g007 
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Figure 8. Influence of corresponding parameters on the pressure and derivative. (A) preferential flow path length L a (B) equivalent 
radius r pa (C) conductivity C da (D) vertical elevation z pa , (E) flow rate constant q DND and (F) storage effect Q> 

doi:1 0.1 371 /journal.pone.0083536.g008 



for a notably long time; therefore, well testing data for this flow 
regime are more likely to be unavailable in the oilfield application 
and not displayed in Fig. 7A. 

The influence of Lp on the pressure transient behavior is shown 
in Fig. 8A. The solutions exhibited in this figure are for a 
preferential flow path with a constant dimensionless conductivity 
of 5 located mid-height in the reservoir [zwD = 0.5), and k = k„is 
also assumed. It indicates that has an effect only upon the early 
stage, and long-time pressure response is unaffected. Fig. 8A 
presents that with the decrease of dimensionless length LD, the 



inertial effect of non-Darcy behavior will cause flow choking and 
increase pressure drop and derivative within the preferential flow 
path. 

The results shown in Fig. 8B indicate that the equivalent radius 
of the preferential flow path rpD has a significant effect on the 
pressure transient performance in the early stage. As the 
dimensionless equivalent radius decreases, the choking effect will 
lead to an increasing pressure drop for a given value of C dD \ the 
smaller the radius, the higher the Pp. For the derivative, as rpi> 
increases, the pressure drop attributable to the choking flow is 
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Table 4. Formation and fluid parameters for Well A1. 



Parameter Name 


Symbol 


Value 


Unit 


Flow rate 


<? 


1 1 1 .84 


m 3 /d 


Fluid density 


P 


850 


kg/m 3 


Fluid viscosity 


f 


0.5 


mPas 


Reservoir porosity 


0 


0.2514 


dimensionless 


Reservoir permeability 


k 


0.2487 


[im 2 


Formation thickness 


h 


18.5 


m 


Compressibility 


Q 


0.001229 


1/MPa 



doi:1 0.1 371 /joumal.pone.0083536.t004 

insignificant compared with the reservoir, indicating a lower value 
of the dimensionless derivative. 

The pressure response shown in Fig. 8C displays the influence of 
dimensionless conductivity. Because the permeability, length and 
radius all have an effect on the conductivity of the preferential flow 
path, the length L D and radius r pD are kept constant. The 
characteristics of transient-pressure behavior for the preferential 
flow path with different conductivities are similar. The effect of 
conductivity concentrates mainly on the early times, yet all the 
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Figure 9. History match of the pressure transient test (A) and 
interewell tracer concentration curve (B) for Well A1. 
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Table 5. Parameters of the preferential flow path for Well A1. 





Parameter Name 


Symbol 


Value 


Unit 


Length of the preferential flow 
path 


L 


46.85 


m 


Permeability of the preferential 
flow path 


k, 


89.42 


|jm 2 


Vertical elevation of the 
preferential flow path 


z p 


7.733 


m 


Equivalent radius of the 
preferential flow path 


<p 


0.7104 


m 


Non-Darcy flow coefficient 


IS 


5.72 x10 6 


m- 1 


Storage constant 


c 


0.024 


m 3 /MPa 



doi:1 0.1 371 /joumal.pone.0083536.t005 

pressure responses coincide at late times. A smaller conductivity 
will contribute to a larger dimensionless pressure and derivative 
when t D -£l. 

Fig. 8D reveals the influence of the dimensionless flow rate 
constant, which represents the non-Darcy effect on the character- 
istics of the pressure transient response. The preferential paths 
with distinct flow rate constants show a similar pressure 
performance, except that the pressure and derivative curves are 
slightly different when t D -SI, As time decreases, the difference 
between these plots increases, which can be attributed to the 
inertial pressure drop caused by non-Darcy flow behavior. The 
comparison with Fig. 8C illustrates that the effect of qnjVD is to 
lower the conductivity of the preferential flow path, which is 
consistent with previous study [70]. The pressure and derivative 
difference for different values of qnjm is detectable at early times 
because the preferential flow path with larger qnjvo behaves as if it 
has a lower conductivity. 

Fig. 8E shows the influence of the vertical location on the 
pressure performance of the preferential flow path. The vertical 
location, represented by the dimensionless distance from the lower 
boundary to the horizontal surface, has a primary effect mainly at 
the early times, whereas the pressure responses at late times 
(formation linear flow) are insensitive to the location. When the 
vertical location is at mid-height in the formation (zpD = 0.5), the 
dimensionless pressure and derivative are smaller. When the 
vertical location is at the bottom (zpD = 0) or top (Zpo = 1) of the 
formation, the pressure and derivative reach the maximum. 
However, if the vertical location is at 1 /4 or 3/4 of the height of 
the formation (zpn =0.25, 0.75), the dimensionless pressure and 
derivative are minimum when tn <0.005. In addition, these 
preferential flow paths located symmetrically in the vertical 



Table 6. Formation and fluid parameters for Well B1. 



Parameter Name 


Symbol 


Value 


Unit 


Flow rate 


q 


30 


m 3 /d 


Fluid density 


p 


798 


kg/m 3 


Fluid viscosity 


/' 


0.44 


mPa-s 


Reservoir porosity 


0 


0.23 


dimensionless 


Reservoir permeability 


k 


0.01609 


[im 2 


Formation thickness 


h 


64.4 


m 


Compressibility 


c, 


0.001173 


1/MPa 



doi:1 0.1 371 /joumal.pone.0083536.t006 
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Figure 10. History match of the pressure transient test (A), 
interewell tracer concentration curve (B) and liquid production 
profile (C) for Well B1. 
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Table 7. Parameters of the preferential flow path for Well B1. 





Parameter Name 


Symbol 


Value 


Unit 


Length of the preferential flow 
path 


L 


56.16 


m 


Permeability of the preferential 
flow path 


k f 


1.394 


jim 2 


Vertical elevation of the 
preferential flow path 


z p 


16.35 


m 


Equivalent radius of the 
preferential flow path 


r P 


0.7568 


m 


Non-Darcy flow coefficient 




2.51 x10 9 


m- 1 


Storage constant 


C 


0.0703 


m 3 /MPa 



doi:1 0.1 371 /journal.pone.0083536.t007 

direction have identical pressure transient behavior, such as Zpn 
= 0.25 and 0.75, z pD =0.1 and 0.9. 

Fig. 8F provides information on the influence of the storage 
effect. The storage affects only the early-time flow period. As Cd 
increases, the dimensionless pressure decreases, and the pure 
storage period is lengthened. For the derivative curve, C D has a 
great influence on the transition period followed by the formation 
linear flow. If the storage is high, the dimensionless pressure 
derivative will increase and display almost flat lines for 0.002 £ t D 
<0.02. 

Characterization of the Preferential Flow Path 

For the purpose of dealing with the preferential flow path by 
successful project design and implementation of conformance 
control treatment, the prerequisites for further Enhanced Oil 
Recovery measures, the characterization of the preferential flow 
path is essential and imperative. As mentioned in the introduction, 
to the knowledge of the authors, non-Darcy flow behavior is 
seldom taken into account in the existing works for the description 
of the preferential path. As shown in Fig. 8 (A-F), all of these type 
curves have similar shapes; therefore, unique solutions are difficult 
to obtain manually. To bridge this gap, an automatic matching 
method has been employed to estimate parameters of the 
preferential flow path because the computer-aided algorithm 
offers numerous advantages over the conventional graphical 
techniques. The automatic history matching technology has such 
merits as the capability to analyze multi-rate tests and actual data 
in the transition regions, to avoid inconsistent interpretations and 
to provide confidence estimates [76]. We use the Gravitational 
Search Algorithm (GSA) [77] to fit the model and the data as 
closely as possible by minimizing the sum of the squares of the 
differences between the measured and model pressures, which is 
given by 



^ ^ \p measured (tk) 
k=\ 



-p mo dei(tk,Cn,kf,li,L,Zp,r p )] (48) 



Conventional gradient-based methods are not efficient and 
practical for solving optimization problems with a high-dimen- 
sional search space because the feasible domain increases 
exponentially with the problem size [77]. If the objective function 
is multimodal and/ or non-smooth, these methods may readily be 
trapped into local optima [78]. To overcome these drawbacks, 
various heuristic approaches, which mimic physical or biological 
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processes in natural phenomena, have been proposed and widely 
adopted in many different areas [79-83]. Among these compu- 
tational intelligence-based techniques, four well-known global 
optimization approaches (Genetic Algorithm (GA) [84], Simulated 
Annealing (SA) [85], Ant Colony Optimization (ACO) [86], and 
Particle Swarm Optimization (PSO) [87]) have attracted the most 
attention and application for their novelty and strong searching 
capacity [88]. However, because of their premature convergence 
and local stagnation, which are frequently observed in many 
applications, these algorithms can achieve a better solution for 
some specific problems rather than others. Searching for novel and 
efficient optimization methods has therefore been an open 
problem [77]. 

Inspired by Newton's law of universal gravitation and mass 
interaction, a new heuristic optimization approach, Gravitational 
Search Algorithm (GSA), was proposed by Rashedi et al. in 2009 
[77]. According to the theory of Newton, each particle in the 
universe attracts every other particle with gravitational force 
directly proportional to the product of their masses [89] . In this 
algorithm, a set of agents, also called masses, are introduced and 
considered as objects to find the optimum. The position of each 
agent denotes a possible solution of the problem, with its mass 
characterizing the fitness value. All of these agents cooperate with 
each other using a direct form of communication through the 
gravitational force, which results in a global movement towards 
other objects with heavier masses. Heavier agents, which 
correspond to better solutions, move more slowly than the lighter 
ones. Ultimately, masses are attracted by the heaviest agents, 
which provide the optimum solution in the search space [88]. 

With a flexible and well-balanced mechanism to enhance both 
exploration (the ability to expand the search space) and 
exploitation (the ability to find the optimum in the vicinity of a 
good solution) abilities, this approach is easy to implement and 
computationally efficient. The high quality performance has been 
certified to solve a wide range of practical optimization problems, 
such as filter modeling [88], future oil demand prediction [90], 
parameters identification [91] and power flow optimization [92]. 
The capacity of Gravitational Search Algorithm has been 
evaluated and compared with the Real Genetic Algorithm 
(RGA), Particle Swarm Optimization (PSO), and Central Force 
Optimization (CFO) for 23 standard benchmark functions, 
including unimodal high-dimensional, multimodal high-dimen- 
sional, and multimodal low-dimensional types. The results show 
that in most cases, GSA provides superior results and in all cases, 
results from GSA are comparable to RGA, PSO and CFO [77]. 
The same conclusion has also been obtained by other researchers 
[93-94]. Therefore, GSA is employed to solve the optimization 
problem and to estimate parameters of the preferential flow path. 
The computational steps of this algorithm can be found in 
[18,77,78,88,91-94]. 

Oilfield Applications 

To demonstrate the validity of the approach described in the 
previous sections, two examples from different oilfields are 
presented. 

Field Case 1: Shengli Oilfield, the third largest reservoir in 
China, is located at the Yellow River Delta in Shandong Province. 
This poorly consolidated sandstone reservoir has been exploited 
for more than 30 years. At present, the field water cut is 96.7%, 
and 34.5% of the oil reserves have been extracted. However, the 
production performance of many injectors/producers and inter- 
well tracer tests reveal that a few preferential flow paths may have 
evolved in the formation. Table 4 lists the formation and fluid 
parameters for Well Al. The observed pressure and derivative are 



plotted in log-log coordinates in Fig. 9A. Using the automatic 
history matching technique of the Gravitational Search Algorithm, 
the computer can obtain the best fit match to the observed data 
with the presented model. Fig. 9A shows an acceptable agreement 
has been discovered, and the interpretation parameters of the 
preferential flow path are given in Table 5. The estimated 
permeability of the preferential flow path is approximately 89.42 
um 2 and of extremely high conductivity, so that the strong 
heterogeneity leads to an incredibly quick premature break- 
through of chemical agents, as observed in the oilfield. 

Polyacrylamide was first injected into the neighbor Well A2, 
300 m from Well Al, to block the high-permeability streak and to 
modify the injection profile. However, only 8.6 hours later, 
polyacrylamide was observed in the production liquid of Well Al. 
With the objective of better understanding the interwell connec- 
tivity, ,9 Fe tracer was injected from the Well A2. Only 7.5 hours 
later, the 59 Fe reached Well Al, and the concentration of 59 Fe 
began to increase significantly. Fig. 9B shows the tracer 
concentration curve for Well Al, and the tracer advance velocity 
is about 40 m/h, which is very high for this conventional reservoir. 
Moreover, integration analysis of geological information, Produc- 
tion Logging Test (PLT) and interwell dynamic connectivity from 
the Capacitance-Resistance Model (CRM) [12] is also consistent 
with the characterization. The treatment project design based on 
this interpretation result also shows very good performance, which 
has demonstrated the validity of this approach. 

Field Case 2: Well Bl belongs to the Huabei Oilfield, located 
in the central area of Hebei Province in China. The formation and 
fluid parameters are shown in Table 6. Automatic history 
matching was also employed to fit the observed pressure and 
derivative with simulated data obtained from our model (Fig. 10A). 
Table 7 lists the estimated parameters of the preferential flow path, 
showing that although the permeability of the ordinary formation 
is only 0.0 1609um 2 , the preferential flow path with permeability as 
high as 1.394 um has evolved around this well. The interwell 
tracer test (Fig. 10B) and liquid production profile (Fig. 10C) also 
verify this conclusion. Considering the reservoir property, the 
tracer advance speed of 2.98 m/d is too high. Fig. 10C also reveals 
that although the net pay of this layer is 64.4 m, 99.1% of the 
production liquid should be attributed to the short segment with a 
thickness of only 1.6 m, which is very close to the estimated 
diameter of the preferential flow path (1.5136 m). 

Conclusions 

Using the revised Forchheimer number recommended by Zeng 
and Grigg [23], the fluid flow behavior within a preferential path 
was examined. The results indicate that non-Darcy flow through a 
preferential path is possible for an actual oilfield and should be 
taken into account for its accurate characterization. This work 
establishes a hybrid analytical/ numerical model for analyzing the 
pressure transient response of the preferential path with non- 
Darcy flow. The characteristics of the pressure transient behavior 
have also been discussed, as well as the influence of each 
corresponding parameter. Based on this model and the Gravita- 
tional Search Algorithm, an automatic matching technology has 
been supplied to estimate parameters of the preferential path with 
non-Darcy flow using the pressure transient test. Oilfield 
applications demonstrate that the characterization is in good 
agreement with interwell surveillance and dynamic analysis. Two 
examples from different oilfields are also provided to substantiate 
the validity of this approach. However, to characterize the 
preferential flow path more accurately and effectively, integration 
analysis with geological information, production logging test and 
several other methods is also essential and recommended. 



PLOS ONE | www.plosone.org 



14 



December 2013 | Volume 8 | Issue 12 | e83536 



Modeling of Preferential Path with Non-Darcy Flow 



Author Contributions 

Analyzed the data: SW Q_F XH. Wrote the paper: SW XH. Conducted the 
theoretical derivation: SW QF. Designed and performed the numerical 
calculation: SW QF XH. 

References 

1. Fattahpour V, Moosavi M, Mchranpour M (2012) An experimental investiga- 
tion on the effect of rock strength and perforation size on sand production. 
J Petrol Sci Eng 86-87: 172-189. 

2. Szymczak P, Ladd AJC (2009) Wormhole formation in dissolving fractures. 
J Geophys Res 114 (B6): B06203. 

3. Colon CFJ, Oelkers EH, SchottJ (2004) Experimental investigation of the effect 
of dissolution on sandstone permeability, porosity and reactive surface area. 
Geochim Cosmochim Ac 68(4): 805-817. 

4. Liu Y, Bai B, Wang Y (20 1 0) Applied technologies and prospects of conformance 
control treatments in China. Oil Gas Sci Tcchnol 65(6): 859—878. 

5. Smith UH, Brigham WE (1965) Field evaluation of watcrflood tracers in a five 
spot. API 

6. Brigham WE, Abbaszadch-Dchghani M (1987) Tracer testing for reservoir 
description. J Petrol Tcchnol 39(5): 519-527. 

7. Beve D, Morrison HF (1991) Borehole-to-surface electrical resistivity monitoring 
of a salt water injection experiment. Geophysics 56(6): 769-777. 

8. Datta-Gupta A, Lake LW, Pope GA, King MJ (1995) A type-curve approach to 
analyzing two-well tracer tests. SPE Form Eval 10(1): 40—48. 

9. Bane RK, Parker RA, Storbeck WG, Sunde RL (1994) Reservoir management 
of the fullerton clearfork unit. In: SPE 27640, presented at Permian Basin Oil 
and Gas Recovery Conference. Midland, Texas, USA; March 16-18. 

10. Li B, Najeh H, Lantz J, Rampurawala MA, Gok I, et al. (2008) Detecting thief 
zones in carbonate reservoirs by integrating borehole images with dynamic 
measurements. In: SPE 1 16286, presented at Annual Technical Conference and 
Exhibition. Denver, Colorado, USA; September 21-24. 

1 1 . Albertoni A, Lake LW (2003) Inferring interwell connectivity only from well-rate 
fluctuations in waterfloods. SPE Reserv Eval Eng 6(1): 6-16. 

12. Sayarpour M, Zuluaga E, Kabir CS, Lake LW (2009) The use of capacitance- 
resistance models for rapid estimation of watcrflood performance and 
optimization. J Petrol Sci Eng 69(3^): 227-238. 

13. Wang S, Jiang H (2010) Determine level of thief zone using fuzzy ISODATA 
clustering method. Transport Porous Med 86(2): 483—490. 

14. Kaviani D, Jensen JL, Lake LW (2012) Estimation of interwell connectivity in 
the case of unmeasured fluctuating bottomhole pressures. J Petrol Sci Eng 90- 
91: 79-95. 

15. Dingcs DD, Ogbc DO (1988) A method for analyzing pulse tests considering 
wellbore storage and skin effects. SPE Form Eval 3(4): 743-750. 

16. Feng Q> Wang S, Gao G, Li C (2010) A new approach to thief zone 
identification based on interference test. J Petrol Sci Eng 75(1-2): 13-18. 

17. Feng Q, Wang S, Wang S, Wang Y, Chen D (201 1) Identification of thief zones 
by dimensionless pressure index in waterfloods. In: SPE 143926, presented at 
SPE Enhanced Oil Recovery Conference. Kuala Lumpur. Malaysia; July 19 21. 

18. Feng Q, Wang S, Zhang W, Song Y, Song S (2013) Characterization of high- 
permeability streak in mature waterflooding reservoirs using pressure transient 
analysis. J Petrol Sci Eng 110: 55-65. 

19. Cui C, Yang Y, Cao G, Zhao X (2009) Effect of high-velocity non-Darcy flow on 
productivity in wormhole of unconsolidated sandstone reservoirs. J Oil Gas 
Tcchnol 31(3): 122-125. 

20. Liu R, Jiang H, Chen M, Dong L (2010) Experiment in high-speed non-Darcy 
flow laws in cross flow channels. Petrol Geol Oilfield Dev Daqing 29(1): 65-69. 

21. Chen Y, Jiang H, Li S (1994) Application of well-to-well tracer test on reservoir 
heterogeneity description. J Univ Petrol 18(supplement): 1-6. 

22. Wang X, Wang J, Wang C, Zeng L, Liu X (2010) Quantitative description of 
characteristics of high-capacity channels in unconsolidated sandstone reservoirs 
using in situ production data. Petrol Sci 7(1): 106-1 1 1. 

23. Zeng Z, Grigg R (2006) A criterion for non-Darcy flow in porous media. 
Transport Porous Med 63(1): 57-69. 

24. Jones SC. (1987) Using the inertial coefficient, p\ to characterize heterogeneity in 
reservoir rock. In: SPE 16949, presented at SPE Annual Technical Conference 
and Exhibition. Dallas, Texas, USA; September 27-30. 

25. Andrade JS, Costa UMS, Almeida MP, Makse HA, Stanley HE (1999) Inertial 
effects on fluid flow through disordered porous media. Phys Rev Lett 82(26): 
5249-5252. 

26. Stanley HE, Andrade JS (2001) Physics of the cigarette filter: tluid How through 
structures with randomly- placed obstacles. Physica A 295(1-2): 17-30. 

27. Forchheimer P (1901) Wasserbewegung durch boden. Zeit Ver Deut Ing 45: 
1782-1788. 

28. Zhang J, Xing H (2012) Numerical modeling of non-Darcy How in near-well 
region of a geothermal reservoir. Gcothcrmics 42: 78—86. 

29. Mahdiyar H, Jamiolahmady M (2011) Improved Darcy and non-Darcy flow 
formulations around hydraulically fractured wells. J Petrol Sci Eng 78(1): 149— 
159. 

30. Wu Y (2002) An approximate analytical solution for non-Darcy flow toward a 
well in fractured media. Water Resour Res 38(3): 1-7. 



31. Perera MSA, Ranjith PG, Choi SK, Airey D (201 1) Numerical simulation of gas 
flow through porous sandstone and its experimental validation. Fuel 90: 547— 
554. 

32. Ghahri P, Jamiolahmady M (2012) A new accurate and simple model for 
calculation of productivity of deviated and highly deviated well-Part I: Single- 
phase incompressible and compressible fluid. Fuel 97: 24—37. 

33. Fourar M, Radilla G, Lenormand K, Moyne C (2004) On the non-linear 
behavior of a laminar single-phase flow through two and three-dimensional 
porous media. Adv Water Resour 27(6): 669-677. 

34. Chai Z, Shi B, LuJ, Guo Z (2010) Non-Darcy flow in disordered porous media: 
A lattice Boltzmann study. Comput Fluids 39(10): 2069-2077. 

35. Chilton TH, Colburn AP (1931) Pressure drop in packed tubes. Ind Eng Chem 
23(8): 913-919. 

36. Fancher G, LewisJA (1933) Flow of simple fluids through porous materials. Ind 
Eng Chem 25(10): 1139-1147. 

37. Tek MR (1957) Development of a generalized Darcy equation. J Petrol Tcchnol 
9(6): 45-47. 

38. Bear J (1972) Dynamics of fluids in porous media. New York: Elsevier. 

39. Scheidegger AE (1974) The physics of flow through porous media (3rd ed.). 
Toronto: University of Toronto Press. 

40. Dybbs A, Edwards RV (1984) A new look at porous media fluid mechanics- 
Darcy to turbulent. In: Bear J, Corapcioglu MY, editors. Fundamentals of 
Transport Phenomena in Porous Media. Dordrecht: Martinus Nijhoff Publisher, 
pp. 199-256. 

41. Blick EF, Civan F (1988) Porous media momentum equation for highly 
accelerated flow. SPE Reserv Eng 3(3): 1048-1052. 

42. Du Plessis JP, Masliyah JH (1988) Mathematical modeling of flow through 
consolidated isotropic porous media. Transport Porous Med 3(2): 145-161. 

43. Green L, Duwez P (1951) Fluid flow through porous metals. J Appl Mech 3(2): 
39-45. 

44. Ma H, Ruth DW (1993) The microscopic analysis of high Forchheimer number 
flow in porous media. Transport Porous Med 13(2): 139-160. 

45. Chukwudozie CP, Tyagi M, Scars SO, White CD (2012) Prediction of non- 
Darcy coefficients for inertial flows through the castlegatc sandstone using 
image-based modeling. Transport Porous Med 95(3): 563—580. 

46. Ergun S (1952) Fluid flow through packed columns. Chem Eng Prog 48(2): 89— 
94. 

47. Thauvin F, Mohanty KK (1998) Network modeling of non-Darcy flow through 
porous media. Transport Porous Med 31(1): 19-37. 

48. ComitiJ, Sabiri NE, Montillet A (2000) Experimental characterization of flow 
regimes in various porous media- III: limit of Darcy 's or creeping flow regime for 
Newtonian and purely viscous non-Newtonian fluids. Chem Eng Sci 55(15): 
3057-3061. 

49. Newman MS, Yin X (201 1) The effect of pore heterogeneity on non-Darcy flow 
by lattice Boltzmann simulation. In: SPE 146689, presented at SPE Annual 
Technical Conference and Exhibition. Denver, Colorado, USA; October 30- 
November 2- 

50. Garrouch AA, Ali L (2001) Predicting the onset of inertial effects in sandstone 
rocks. Transport Porous Med 44(3): 487-505. 

51. Belhaj HA, Agha KR, Nouri AM, Butt SD, Islam MR (2003) Numerical and 
experimental modeling of non-Darcy flow in porous media. In: SPE 81037, 
presented at SPE Latin American and Caribbean Petroleum Engineering 
Conference. Port-of-Spain; April 27-30. 

52. JanicekJD, Katz DL (1955) Applications of unsteady state gas flow calculations. 
In: Presented at Research Conference on "Flow of Natural Gas from 
Reservoirs". University of Michigan; June 30-July 1. 

53. GeertsmaJ (1974) Estimating the coefficient of inertial resistance in fluid flow 
through porous media. SPE J 14(5): 445-450. 

54. MacDonald IE (1979) Flow through porous mcdia-the Ergun equation revisited. 
Ind Eng Chem Fund 18: 189-208. 

55. Pascal H, Quillian RG, Kingston J (1980) Analysis of vertical fracture length and 
non-Darcy flow coefficient using variable rate tests. In: SPE 9348, presented at 
SPE Annual Technical Conference and Exhibition. Dallas, Texas, USA; 
September 21-24. 

56. Coles ME, Hartman KJ (1998) Non-Darcy measurements in dry core and the 
effect of immobile liquid. In: SPE 39977, presented at SPE Gas Technology 
Symposium. Calgary, Alberta, Canada; March 15-18. 

57. Li D, Svec RK, Engler TW, Grigg RB (2001) Modeling and simulation of the 
wafer non-Darcy flow experiments. In: SPE 68822, presented at SPE Western 
Regional Meeting. Bakcrsficld, California, USA; March 26-30. 

58. Friedel T, Voigt HD (2006) Investigation of non-Darcy flow in tight-gas 
reservoirs with fractured wells. J Petrol Sci Eng 54: 112-128. 

59. Macini P, Mesini E, Viola R (2011) Laboratory measurements of non-Darcy 
flow coefficients in natural and artificial unconsolidated porous media. J Petrol 
Sci Eng 77(3-4): 365-374. 



PLOS ONE | www.plosone.org 



15 



December 2013 | Volume 8 | Issue 12 | e83536 



Modeling of Preferential Path with Non-Darcy Flow 



60. Ahmed T, McKinncy PD (2005) Advanced reservoir engineering. Burlington: 
Gulf Professional Publishing. 

61. Chen C, SongX, Li J (2012) Dominant How channels of point-bar reservoirs and 
their control on the distribution of remaining oils. Acta Petrol Sin 33(2): 257- 
263. 

62. Gringarten AG, Ramey Jr HJ (1974) Unsteady-state pressure distributions 
created by a well with a single horizontal fracture, partial penetration, or 
restricted entry. SPE J 14: 413^126 

63. Ehlig-Economides GA, Joseph J (1987) A new test for determination of 
individual layer properties in a multilayered reservoir. SPE Form Eval 2: 261- 
283. 

64. Kuchuk EJ, Gnur M, Hollaender F (2010) Pressure transient formation and well 
testing: convolution, deconvolution and nonlinear estimation. Amsterdam: 
Elsevier. 

65. Ezulike O, Igbokoyi A (2012) Horizontal well pressure transient analysis in 
anisotropic composite reservoirs-A three dimensional semi- analytical approach. 
J Petrol Sci Eng 96-97: 120-139. 

66. Liu X, Zhao G (2005) Transient -pressure behavior of cold heavy-oil production 
wells. In: SPE 97791, presented at International Thermal Operations and Heavy 
Oil Symposium. Calgary, Alberta, Canada, November 1-3. 

67. Guppy KH, Ginco-Lcy H, Ramey HJ (1981) Effect of non-Darcy How on the 
constant-pressure production of fractured wells. SPEJ 21(3): 390—400. 

68. Guppy KH, Cinco-Ley H, Ramey HJ, Samaniego-V F (1982) Non-Darcy flow 
in wells with finite -conductivity vertical fractures. SPEJ 22(5): 681-698. 

69. Gringarten AC, Ramey HJ (1973) The use of source and Green's functions in 
solving unsteady-flow problems in reservoirs. SPEJ 13(5): 285-296. 

70. Ozkan E, Sarica C, Haciislamogln M, Raghavan R (1995) Effect of conductivity 
on horizontal well pressure behavior. SPE Adv Technol 3(1): 85-94. 

71. Agarwal RG, Al-Hussainy R, Ramey HJ (1970) An investigation of wellbore 
storage and skin effect in unsteady liquid flow: I. Analytical treatment. SPE J 
10(3): 279-290. 

72- Cinco-Ley H, Samaniego-V F (1977) Effect of wellbore storage and damage on 
the transient pressure behavior of vertically fractured wells. In: SPE 6752, 
presented at SPE Annual Technical Conference and Exhibition. Denver, 
Colorado, USA; October 9-12. 

73. Al-Kobaisi M, Ozkan E, Kazemi H (2006) A hybrid numerical/ analytical model 
of a finite -conductivity vertical fracture intercepted by a horizontal well. SPE 
Reserv Eval Eng 9(4): 345-355. 

74. Cinco-Lcy H, Samaniego-V F (1981) Transient pressure analysis for fractured 
wells. J Petrol Technol 33(9): 1749-1766. 

75. Ozkan E, Raghavan R, Joshi SD (1989) Horizontal-well pressure analysis. SPE 
Form Eval 4(4): 567-575. 

76. Home RN (1990) Modern well test analysis: A computer-aided approach. Palo 
Alto: Petroway Inc. 



77. Rashedi E, Nczamabadi-pour H, Saryazdi S (2009) GSA: A gravitational search 
algorithm. Inf Sci 179: 2232-2248. 

78. Shaw B, Mukherjee V, Ghoshal SP (2014) Solution of reactive power dispatch of 
power systems by an opposition-based gravitational search algorithm. Int J Elee 
Power 55: 29-40. 

79. Hajizadeh Y, Christie M, Demyanov V (2011) Ant colony optimization for 
history matching and uncertainty quantification of reservoir models. J Petrol Sci 
Eng 77: 78-92. 

80. Al-Gosayir M, Babadagli T, Leung J (2012) Optimization of SAGD and solvent 
additive SAGD applications: Comparative analysis of optimization techniques 
with improved algorithm configuration. J Petrol Sei Eng 98-99: 61—68. 

81. Li X, Yin M (2012) Application of differential evolution algorithm on self- 
potential data. PLoS ONE 7(12): c51 199. 

82. Gao N, Yang N, Tang J (2013) Ancestral genome inference using a genetic 
algorithm approach. PLoS ONE 8(5): c62156. 

83. Hamadnch N, Khan WA, Sathasivam S, Ong HC (2013) Design optimization of 
pin fin geometry using particle swarm optimization algorithm. PLoS ONE 8(5): 
c66080. 

84. Tang K, Man K, Kwong S, He Q, (1996) Genetic algorithms and their 
applications. IEEE Signal Proc Mag 13(6): 22-37. 

85. Kirkpatrick S, Gclatto CD, Vecchi MP (1983) Optimization by simulated 
annealing. Science 220: 671-680. 

86. Dorigo M, Maniezzo V, Colorni A (1996) The ant system: optimization by a 
colony of cooperating agents. IEEE Trans Syst, Man, and Gybern-Part B 26(1): 
29^1. 

87. Kennedy J, Eberhart RC (1995) Particle swarm optimization. Proceedings of 
IEEE International Conference on Neural Networks 4: 1942-1948. 

88. Rashedi E, Nczamabadi-pour H, Saryazdi S (2011) Filter modeling using 
gravitational search algorithm. Eng Appl Artiflntel 24: 117—122. 

89. Halliday D, Resnick R, Walker J (1993) Fundamentals of physics. New York: 
John Wiley and Sons. 

90. Behrang MA, Assareh E, Ghalambaz M, Assari MR, Noghrehabadi AR (201 1) 
Forecasting future oil demand in Iran using GSA (Gravitational Search 
Algorithm). Energy 36(9): 5649-5654. 

91. Li C, Zhou J (2011) Parameters identification of hydraulic turbine governing 
system using gravitational search algorithm. Energ Gonvers Manage 52: 374— 
381. 

92. Duman S, Guvenc U, Sonmcz Y, Yorukeren N (20 1 2) Optimal power flow using 
gravitational search algorithm. Energ Convers Manage 59: 86-95. 

93. Chatterjee A, Mahanti GK, Pathak N (2010) Comparative performance of 
gravitational search algorithm and modified particle swarm optimization 
algorithm for synthesis of thinned scanned concentric ring array antenna. Prog 
Electromagn Res B 25: 331-348. 

94. Eldos T, Qasim R (2013) On the performance of the Gravitational Search 
Algorithm. Int J Adv Comput Sci Appl 4: 74-78. 



PLOS ONE | www.plosone.org 



16 



December 2013 | Volume 8 | Issue 12 | e83536 



